The Olami-Feder-Christensen earthquake model in one dimension 
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We study the earthquake model by Olami, Feder and Christensen in one dimension. While the 
size distribution of earthquakes resembles a power law for small system sizes, it splits for larger 
system sizes into two parts, one comprising small avalanches and showing a size independent cutoff, 
and the other comprising avalanches of the order of the system size. We identify four different types 
of attractors of the dynamics of the system which already exist for very small systems. For larger 
system sizes, these attractors contain large synchronized regions. 

PACS numbers: 05.65.+b,45.70.Ht 



I. INTRODUCTION 

The Olami-Feder-Christensen earthquake model [l| is 
probably the most studied nonconservative and suppos- 
edly self-organized critical model. Nevertheless, the ori- 
gin of its power-law like avalanche-size distribution is still 
not clear. Apart from these power laws, the model shows 
a variety of other interesting and unusual features such as 
a marginal synchronization of neighbouring sites driven 
by the open boundary conditions 0, and the violation 
of finite-size scaling Q| together with a qualitative 
difference between system- wide earthquakes and smaller 
earthquakes 5j. Also, small changes in the model rules 
(like replacing open boundary conditions with periodic 
boundary conditions 6], or introducing frozen noise Q), 
destroy the power laws. Recently, it was found that the 
results of computer simulations are strongly affected by 
the computing precision § , and that the model exhibits 
sequences of foreshocks and aftershocks . 

In order to better understand the model, we study here 
its one-dimensional version. The model is highly non- 
trivial even in one dimension, and some of its properties 
resemble those in two dimensions. Just as in two dimen- 
sions, we find large synchronized regions and a funda- 
mental difference between the avalanches triggered at the 
boundaries and those triggered deep inside the system. 
We identify different types of attractors of the dynam- 
ics of the system and explain the features of the model 
in terms of the properties of these attractors. Our main 
finding is that the system in the stationary state can 
be separated into a boundary region, where all larger 
avalanches are triggered, and one (or two) synchronized 
inner regions, the size of which can be varied without 
changing the behaviour of the boundary region. 

The outline of this article is as follows. In the next 
section, we introduce the model rules. In section ITTT1 we 
focus on a system of up to 4 lattice sites and find its at- 
tractors. In section lTVl we view the model from a dynam- 
ical systems' perspective and present a general analytical 
approach that allows to classify the attractors into four 
different types. In section [V] we study larger systems. 
First, we investigate the approach to the stationary state 
as function of the system size and the model parameter. 
Then, we discuss the properties of the stationary state. 



Finally, we summarize and discuss our main findings in 
section Ell 



II. THE MODEL 

The Olami-Feder Christensen model is a discretized 
and simplified version of the Burridge-Knopoff model of 
earthquakes |lOj| . In an one-dimensional system consist- 
ing of L sites, it is defined by the following rules: At each 
site i = 1, . . . , L, a continuous variable Zi is defined that 
represents a local force. The force at all sites increases 
uniformly at a constant rate, which we set equal to 1. 
When the force exceeds the threshold value z c , which 
can be chosen to be z c — 1 without loss of generality, the 
force at this site is reset to zero, while the two nearest 
neighbours (or the only neighbour, if the toppling site is 
a boundary site) receive a force increment of azi. The 
parameter a is the only parameter of the model, and it 
has a value in the interval (0, 0.5). If a neighbour is lifted 
above the threshold, the force on its neighbours is imme- 
diately increased according to the same rule, etc., until 
the "avalanche" (the earthquake) is finished. The "size" 
of an avalanche, s, is denned to be the number of top- 
pling events during this avalanche. Such an earthquake 
is instantaneous on the time scale of driving. After the 
earthquake, the force is again increased with unit rate, 
until the next site reaches the threshold, triggering the 
next earthquake, and so on. 

This model is deterministic, once the initial conditions 
are given. Usually, the initial conditions are chosen ran- 
domly from a uniform probability distribution for each 
site. Since the model is deterministic and dissipative, it 
has attractors of the dynamics. 

From a dynamical systems' perspective, the model can 
be viewed as a L-dimensional map, which maps the state 
of the model after one avalanche (which may have size 
1) on the state after the next avalanche. Due to the 
toppling, the map is non continuous. 

If a > a c = (\/3 — 1) /2 ~ 0.366, a site that topples can 
in principle receive from its neighbours packages of a total 
size larger than 1, causing the first site to topple again. 
Throughout this paper, we assume that each site topples 
only once during an avalanche, and we limit therefore our 
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numerical studies and analytical arguments to the case 
a < a c , except where indicated otherwise. 

Let us briefly summarize a few known results that 
are relevant for our study of the one-dimensional sys- 
tem. Firstly, it is found that for periodic boundary con- 
ditions the dynamics approach an attractor where every 
site gives to its neighbours only force packages of size 
a. This means that no force value exceeds the threshold 
value z c on the attractor. After a time 1 — 2a, each site 
has toppled once and has received two force packages of 
size a from its neighbours. This means that after time 
1 — 2a the force on each site is again the same. Slightly 
increasing or decreasing force values gives again a peri- 
odic orbit with period 1 — 2a, as long as this change does 
not cause a toppling site to lift its neighbour above the 
threshold. 

Secondly, the open boundary conditions are responsi- 
ble for the occurrence of large avalanches and large syn- 
chronized regions, where neighbouring sites differ in force 
values only by a small amount. A nice explanation for 
this has been suggested by Middleton and Tang ten years 
a S° 0- They considered a system of two sites where 
one site is driven at a slower rate than the other. This 
mimicks the fact that sites close to the boundary receive 
on an average less force packages than those deep inside 
the system. The two sites settle on an attractor where 
the slower site always topples first and lifts the faster 
site above the threshold. The faster site therefore looses 
more force when it topples, and the slower site receives 
a larger force package. This compensates the different 
driving speed, and the two sites remain synchronized and 
always topple together. 

Thirdly, the largest possible force package that a site 
can pass to its neighbour is a/(l — a). This package size 
is reached if an avalanche passes through a region where 
all sites are at the threshold. 



III. VERY SMALL SYSTEMS 
A. L = 2 

For L = 2, any state where the force difference between 
the two neighbours is larger than a, is part of a cycle of 
period 1 — a. After a time 1 — a, each site has toppled 
once and has received a force package of size a due to the 
toppling of the neighbour. Furthermore, it has received 
a force increment 1 — a due to the uniform driving. The 
net change in force is therefore zero. 



B. L = 3 

For L = 3, we consider the state of the system when- 
ever the middle site has toppled, i.e., when zi = 0. z\ 
and Z3 then take values in [a, 1 + a + a 2 ] only. The lower 
value a is realized, if a site just toppled itself and lifted 



Z2 exactly to the threshold. The upper limit occurs only 
if all sites were at the threshold before. 

Without loss of generality, we assume z\ < Z3. Dif- 
ferent regions in the Z1—Z2 plane can be characterized 
by their sequence of topplings t% and of growth, until 
22 = again. An example of this would be (t$, <7, ti, ^2), 
for which first the right site topples, followed by a uni- 
form growth, and then the toppling of the left site lifts zi 
over the threshold. There is a total of 14 different such 
regions, seven of which are marked by letters a to g in 
figure 

By investigating the transitions between these regions, 
one finds two attractive fixed point lines. One fixed point 
line is at z x = z* = " ( 1 1 +^ a) with z 3 E (z*, 1]. The other 
fixed point line is obtained by interchanging Z\ and Z3. 
The corresponding attractor is a cycle of two avalanches, 
written as (g, tj, g, ti, t^) in the above notation, with i be- 
ing the site with = z* , and j being the other neighbour 
of site 2. 

A special case is the symmetric case z\ = 23 = z*. Site 
1 and 3 topple simultanously, thereby lifting site 2 above 
the threshold, and receiving a package z* in turn. 

To summarize, a system with L = 3 approaches a pe- 
riodic attractor with 2 different avalanches. 
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FIG. 1: Different regions in the 21-23 plane (with 22 = 0) 
that have different toppling sequences. Only the part with 
z\ < 23 is shown. Thin lines separate regions, the thick line 
is the fixed point line. 



C. L = 4 

For the system size L = 4, all attractors of the dynam- 
ics are periodic. We find a variety of different attractors 
for a given value of a. Figure shows the period of the 
attractors found doing a scan of the same 128 4 different 
initial conditions for each value of a and for 2 different 
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precisions. The period is measured once in terms of the 
number of topplings, # t , and once in terms of the number 
of avalanches, # a . 
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FIG. 2: Period of attractors as function of a, measured in 
the number of avalanches per period # a (dots) and in the 
number of topplings per period #t (circles). The solid line 
corresponds to f(x) = 4/x. Top: precision 2 -20 ; bottom: 
precision 2 -10 . The points above 10 6 for a = 0.34 are valid 
data points. 

One can discern the following features: 

1. Degeneracy: Attractors with different numbers of 
avalanches per period have the same total number 
of topplings per period. We found that different 
attractors can have a different toppling sequence, 
while the force packages that each sites receives 
from its neighbours are identical. One explanation 
for this is that there are mechanisms in the model 
that create a degeneracy of different sites, for in- 
stance when two sites have been the end points of 
the same avalanche (after which they both had zero 
force) and remain synchronized until they reach the 
threshold again. The toppling sequence then de- 
pends on the exact implementation of the algorithm 
and on rounding errors due to finite precision; one 



or the other version of the attractor can be ob- 
tained depending on the initial conditions. The to- 
tal number of topplings and the packages sizes can 
be identical on the two versions of the attractor if 
there is a site between the two synchronized sites 
that topples only after having received a package 
from both sides. 

2. Persistence: Attractors exist over a certain inter- 
val of a values. If a is changed slowly (such that 
the system can follow adiabatically) , the attractor 
remains the same as long as the avalanches remain 
the same. Eventually, a point will be reached where 
an avalanche decays into two avalanches (because 
the size of a package is no more large enough to lift 
the neighbour over the threshold), or where two 
avalanches merge to form one avalanche (because 
the distance between two neighbours has become 
smaller than the package size). As we will see in 
the next section, the stability of the attractor can 
change at this point (but typically not before). If 
the state with the new toppling pattern is stable, 
we have a new attractor with the same period (if 
measured in number of topplings). Otherwise, if 
the new state is unstable, the system moves to a 
different attractor, and we obtain a step in Figure 
□ 

3. Divergence of the periods at small a: At small a, 
the period of the attractors is increasing. This can 
be explained by considering a boundary site and 
its neighbour. During each time interval 1 — 2a, 
the boundary site receives a package of the order 
of a from its neighbour, while the neighbour re- 
ceives two packages of the order of a. It therefore 
takes of the order of 1/a time intervals until the 
two sites have again roughly the same height. For 
this reason, the period diverges as 1/a. 

4. Smallest possible period: For all values of a, the 
shortest attractor has 4 topplings. In configuration 
space, this corresponds to the state (z, z, 0, a + a 2 ), 
with an arbitrary force z£ [a + a 2 , 1] , and with the 
toppling sequence g, t\, t 2 , (g), t±, t 3 . Whether this 
attractor is realized, depends on the implementa- 
tion of the algorithm. For smaller a and for smaller 
precision, it occurs more often and eventually has 
the weight one. 

5. Vastly different periods: For a given a, there ex- 
ist attractors with widely different periods. The 
most prominent periods lie in two bands, which are 
clearly visible in the figure. For certain values of a, 
very large periods occur with a considerable weight. 
Attractors with these large periods typically have 
toppling sequences that are most of the time peri- 
odic with a much shorter period than that of the at- 
tractor, but the force values do not have this short 
period. 
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Sensitivity to the computing precision: Attractors 
with larger periods occur (for the same value of a) 
less often when the computing precision is smaller. 
The reason for this is that on longer attractors, 
there are more states at all and states that are close 
to each other are more likely to occur. In a simu- 
lation with smaller numerical precision, such states 
can become identical, and the period of the attrac- 
tor becomes shorter. 



IV. ANALYTICAL APPROACH AND 
CLASSIFICATION OF ATTRACTORS 

We describe the state of the system as a difference vec- 
tor x = (x\, . . . , xl-x) with Xi — Zi+i — Zi. The uniform 
increase in force does not change this vector, but top- 
pling sites do change it. The force value of the first site 
toppling in an avalanche is decreased by one, while its 
two neighbours receive a package a. This process can be 
decribed by adding a vector 



(0, . . . , 0, —a, 1 + a, —(1 + a), a, 0, 



• 0) 



to x (with the four nonzero elements at the appropriate 
place). In contrast, a subsequent toppling event can be 
described by applying to 5r a matrix that is identical 
to the unity matrix everywhere except for a 4 x 4 block 
on the diagonal. There are two different matrices, corre- 
sponding to avalanches propagating to the right and to 
the left. We assume that site i just toppled. Then, the 
next toppling event is represented by the matrix 
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(1) 



if the avalanche is moving to the left. The nontrivial 
column (number i— 1) describes the toppling of site (i— 1), 
which was lifted above the threshold by the prior toppling 
of site i. 

If the avalanche is moving to the right, the correspond- 
ing matrix reads 



/I 



Mi +1 = 



1 a 

-a 

1 + a 1 

-a 1 







(2) 



with the nontrivial column being number i. If site i — 1 
and site i + 1 both topple simultaneously, the matrix 
contains both nontrivial columns. If the toppling site is 
a boundary site or the site next to it, surplus columns 
and rows of the nontrivial block of or M* arc to be 
removed. 

We now focus our interest on the difference 5x between 
two systems and assume that the two systems have the 
same toppling sequence. This will be the case as long as 
they are sufficiently close to each other. Then they will 
both be updated by adding the same vectors and multi- 
plying the same matrices in the same order. Adding the 
same vector to the state x of both systems has no effect on 
the difference Sx. The difference between the states of the 
two systems will therefore evolve solely by multiplying 
matrices of the form MJ and M* to it (- except if the first 
toppling site lifts both neighbours above the threshold; in 
this case the first matrix associated with this avalanche 
contains two nontrivial columns, while the other matri- 
ces have the usual form, since the two branches of the 
avalanche commute after the first toppling.) 

Let one of the two systems be on a periodic orbit. 
Whether the other system will approach the orbit, de- 
pends on the largest eigenvalue of the product 



s = 

v(P) v{P- 



1) 



»(2) »*»(!) 



.MT'M 



K2) iTi Ki) 
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with P being the total number of matrices occurring dur- 
ing one period, i(p) being the nontrivial column index 
for matrix number p, and v{p) being I or r depending 
on whether an avalanche is moving to the left or to the 
right. 

If the largest eigenvalue of S is larger than one, the or- 
bit is unstable and cannot be an attractor of the dynam- 
ics. In the following, we describe four types of attractors 
that we have found in the model. All these attractors 
occur already in a small system of L = 4, but are also 
seen in large systems. 



A. Marginally stable attractors 

Marginally stable attractors occur if the largest eigen- 
value of S is identical to 1. We observed regularly attrac- 
tors where a column i of S is identical to the unit vector 
e*i . This means that e*j is an eigenvector of S , and adding 
a small multiple ee* to the periodic orbit gives again a pe- 
riodic orbit with the same toppling sequence. In terms of 
forces Zj this means that increasing or decreasing all force 
values Zj with j < i by a small amount results again in a 
periodic orbit. The product 10 contains no matrix MJ or 
M* +1 . Sites i + 1 and i never cause each other to topple. 
There is no avalanche that includes simultaneously site 
i + 1 and site i. We say that there is a "barrier" between 
sites i + 1 and i. We found that the total size of the force 
packages that site i + 1 gives to site i during one period is 
identical to the total size of the force packages that site 
i gives to site i + 1. For L = 4, the barrier is always in 
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the middle of the system. For larger sizes, it need not 
be in the middle, but it is often found at the center of 
a system, since the synchronization proceeds at constant 
speed from the boundaries (see below). 

Sites i and i + 1 to the right and left of the barrier 
must topple equally often during one period. If this was 
not the case and if site i toppled more often, there would 
be an instance where site i topples twice without site 
z + l toppling in between. After site i has toppled for 
the first time, its force is zero, and that of site i + 1 is 
at least as large as a. In order for site i to reach the 
threshold before site i + l, it must receive a package from 
its left neighbour that is larger than a, while site i + l 
receives no package. The largest possible package size 
is a/(l — a), and therefore site i + l has at least the 
force 1 — a 2 /(l — a) > 1 — a at the moment where site 
i reaches the threshold for the second time. This means 
that site i+ 1 is lifted above the threshold by the toppling 
of site i, in contradiction to our assumption that there is 
a barrier between the two sites. Therefore, the two sites 
must topple equally often. 

There exists no attractor with 2 or more barriers. We 
show this in two steps. First, let us assume that all force 
packages passed over the barriers are of size a. Then 
the region between the two barriers is like a system with 
periodic boundary conditions, and no package passed on 
within this region is larger than a. The two sites imme- 
diately outside the barriers must not be lifted above the 
threshold by their neighbours. Otherwise, they would 
pass packages larger than a over the barrier. Further- 
more, the two sites immediately outside the barriers must 
not topple more often than the sites between the barri- 
ers. Therefore they cannot receive packages larger than a 
from their outward neighbours. This means that the sites 
immediately outside the barriers in fact also belong to the 
domain of sites that are never lifted above the threshold 
and that always receive packages of size a. By repeating 
this argument, we find that no site in the system can be 
lifted above the threshold. However, this situation can- 
not be realized with open boundary conditions. It occurs 
for periodic boundary conditions. 

Now, since we have ruled out the possibility that the 
region between the two barriers receives only force pack- 
ages of size a from outside, let us assume next that they 
receive on an average packages larger than a. We sim- 
ulated the region between the two barriers by inserting 
packages of a size larger than a at its boundaries imme- 
diately after the boundary sites have toppled. This leads 
to attractors where avalanches are triggered at the center 
of the system and are running outwards. The attractors 
and therefore the number of topplings per unit time of 
the boundary site are determined by the size of the re- 
gion and the size of the packages received from outside. 
On the other hand, this number of topplings must be 
identical to the number of topplings of the site on the 
other side of the barrier in the original system. However, 
there is no free continuous parameter left to match this 
condition, and it can therefore usually not be satisfied. 



This problem does not arise in the case of a single barrier, 
because the state of the system can be symmetric about 
the barrier, thus satisfying the matching conditions. 

Finally, let us consider a periodic orbit at the bound- 
ary of the basin of attraction of the marginally stable 
attractors. In order to obtain this orbit, we increase or 
decrease all force values Zj with j < i by an amount 
such that there is a moment in time where site i (or 
i + l) is lifted by site i + l (or i) exactly to the thresh- 
old. The metastable orbit has now become degenerate 
with an orbit where site i (or i + 1) is lifted by site i + l 
(or i) infinitcsimally above the threshold. This orbit has 
no barrier, and the matrix S corresponding to this peri- 
odic orbit is different from the one corresponding to the 
metastable orbit. Its largest eigenvalue will therefore be 
different from 1. We have seen realizations of the inter- 
esting case that the largest eigenvalue becomes smaller 
than 1. This means that the periodic orbit that is at 
the boundary of the basin of attraction of the marginally 
stable attractors can itself be an attractor that is reached 
from a nonvanishing set of initial conditions. 



B. Strongly stable attractors 

In general, the matrix S is obtained by starting with 
the unit matrix and multiplying the appropropriatc ma- 
trices M; or one after the other. 

The effect of (or ) on a matrix is that row num- 
ber i (or row number i — 1), multiplied by a certain factor 
±a or (1 + a), is added to 3 neighbouring rows and is 
itself multiplied by —a. 

Let us perform an expansion in powers of a and focus 
on the elements of order a . In order a , the matrices 
M| and M* simply add one row to a neighbouring one 
and replace the original row with all zeros. If a boundary 
site is caused to topple by its neighbour, all l's that have 
been in the boundary row are flushed out of the system. 
Starting with a unit matrix, all l's that remain in the 
system are in the same row after multiplying a sufficient 
number of M.\, matrices. 

We have never seen an attractor where this does not 
happen. Any given site of the system is reached by an 
avalanche that starts near the boundary, and therefore 
there cannot be l's left in different rows. The marginally 
stable attractors have all l's in a row i, where they stay 
forever. 

However, there is also the possibility that all the l's 
arc flushed out of the system by avalanches that extend 
from inside the system to the boundary. In this case, the 
largest eigenvalue of S is of the order a, and the attractor 
is quickly approached. We have seen many examples of 
such strongly stable attractors. 
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C. Weakly stable attractors 

If not all l's are flushed out of the system, and if there 
is no barrier in the system, the largest eigenvalue of S 
belonging to an attractor is 1 — 0(a n ) with some power 
n of a. This corresponds to the situation where the row 
containing the l's remains in the system and is moved 
around by avalanches coming from both directions. 

If a is small or n is large, attractors are approached 
very slowly. For a — 0.2, we have seen attractors with 
n = 2 for L = 4 and an estimated n = 11 for L = 
20. However, we did not attempt a systematic survey of 
relaxation times towards the attractors as function of a 
and L. 
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D. Complex attractors 

There exist attractors with strikingly long periods of 
the order of many thousands for L = 4 or much larger 
for larger L (if period is measured in total force incre- 
ment per site). Typically, these attractors contain long 
quasiperiodic sections where the sequence of avalanches 
remains the same but the force values change slowly, just 
as one would expect close to a weakly stable or weakly 
unstable periodic orbit. This quasiperiodic sequence is 
eventually interrupted by an intermittent phase contain- 
ing other avalanches, until the quasiperiodic phase is en- 
tered again. One can understand the origin of such com- 
plex attractors in the following way. Imagine a weakly 
stable attractor for a certain value of a. Now change 
a slowly and let the system follow adiabatically. The 
largest eigenvalue of S on the resulting attractor will have 
the same coefficients if expanded in powers of a, as long 
as the avalanches remain the same. Eventually, a value 
of a will be reached where two avalanches merge or an 
avalanche splits, changing the product of M matrices, 
which now can have an eigenvalue larger than 1. Never- 
theless, there may still be a region nearby in state space 
where the old sequence of avalanches can be maintained 
for a long time if the largest eigenvalue of S of the old 
attractor was only slightly smaller than 1 (i.e. if the old 
attractor was weakly stable). 

V. NUMERICAL INVESTIGATION OF LARGE 
SYSTEMS 

A. Transient stage 

We now present simulation results for larger systems. 
We first study the transition from a random initial state 
to a stationary state. Figure shows the force values 
throughout a system of size L = 1000 for different times. 
One can see that starting from the boundary more and 
more sites become synchronized, until in the stationary 
state all sites apart from a few sites at the boundary have 
almost the same force z,. 



FIG. 3: Random initial configuration and after 3 x 10 6 , 5 x 10 6 
and 9 x 10 6 avalanches in configuration space for L = 1000 
and q = 0.2 (from top left to bottom right). 

We have evaluated the transient time using two differ- 
ent measures of the degree of synchronization: (i) The 
standard deviation 

a 2 (t)=mm t (z 2 (r)-^(r)) 

as function of a and L, where the bar denotes the average 
taken over all sites i. (ii) Nearest-neighbour deviation 

crjvjv(i) = min (znn{t) - z 2 (t)) 

0<T<t \ / 

where again the bar denotes the average over all i and 
znn = jfzj Ej=i z i z i+i Time was measured as the 
number of topplings per site (i.e. the total number of 
topplings divided by the system size). 

Figure ^ shows our results for these two synchroniza- 
tion measures, averaged over 1000 random initial config- 
urations. The two figures at the top show the behaviour 
of <jnn (left) and a (right) as a function of time and of 
the system size L. The lower plateau indicates clearly 
the stationary state. One finds that the transient time 
is proportional to the system size. This is due to the in- 
ward proceeding synchronization, which takes place at a 
constant rate, if L is sufficiently large. For small L, the 
boundary layer takes a large part of the system, and there 
is therefore little synchronization. Apart from the tran- 
sition to the stationary state with a small value of a and 
&nn, one can also distinguish an earlier transition, where 
the two measures leave the value 1/12 corresponding to 
a random initial configuration. We interpret this transi- 
tion as the onset of the formation of synchronized blocks, 
after the boundary layer has been set up. The character- 
istic shape of the curves between these two transitions is 
strikingly different. While the nearest-neighbour devia- 
tion decreases rather fast once the synchronization starts, 
the standard deviation remains on a second plateau until 
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FIG. 4: Standard deviation a (right) and nearest neighbour 
deviation ctjvjv (left) for fixed value of a — 0.35 (top) and fixed 
value of L — 200 (bottom) as function of time (measured in 
number of topplings per site) and L and a respectively, note 
that time is plotted in log-scale and the double logarithmic 
scale in the upper two figures. 



its final decrease. This comes from the two synchronized 
blocks, which usually have a different force value, as one 
can see in the second picture of figure [3] 

A similar behaviour is found for the dependency on a 
as shown in the lower two figures for a fixed system size 
L = 200. For values of a below the critical value a c , the 
nearest-neighbour deviation depends only weakly on a. 
The sharp transition at the end of the high plateau of 
(Jnn is linear in a. The bulky part above a c must be 
due to the fact that a site can now topple twice during 
the same avalanche. While the onset of synchronization 
is better visible in the data for <jnn, which decay very 
rapidly, the transition to the stationary state is much 
better visible in the data for a, which remain close to 
the initial value for a longer time. A more detailed in- 
vestigation of the transition time to the stationary state 
reveals the following: (i) Over a wide range of a values 
this transition time depends exponentially on a, as shown 
in Figure El where the time is plotted that is needed to 
reach a — 0.01. (ii) For small values of a, the data show 
a power law with an exponent around —2.84 (Figure El . 
Since the synchronization proceeds very slowly, we did 
not measure the time to the stationary state, but the 
time to show for the first time an avalanche larger than 
20. The following analytical arguments suggest that 
the transient time should indeed diverge at least as fast 
as aT 2 with a. Let us define a time unit as the time 
during which a force 1 — 2a is added to the system. A 
site at the center receives two packages of size a from its 
neighbours per unit time and topples on an average once 
per unit time. A boundary site receives only one package 
and has on an average approximately 1 — a topplings per 
unit time. A site in the synchronized block topples on an 
average y times per unit time, with y being intermediate 




0.1 0.2 0.3 0.4 0.5 

a 



FIG. 5: Time (measured in topplings per site) needed to 
reach a = 0.01 for a system of L — 200 sites as a function of 
a, averaged over 1000 different random initial systems. The 
peaks are real and are not statistical fluctuations. 
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FIG. 6: Time (measured in topplings per site) needed to 
generate the first avalanche of size > 21 as function of a, av- 
eraged over 1000 different random initial states. The straight 
line corresponds to /(a) oc a~ 2 ' 84 . 

between these two limit cases, 1 — a < y < 1. Initially, 
the force difference between the synchronized block and 
the site that will be synchronized next, is of the order 
of 1. In order to decrease this force difference to a value 
of the order of a, the difference in the total number of 
topplings between the block and its neighbour must be 
of the order of l/a, which is achieved after a time of the 
order of l/a(l — y). This increases with decreasing a at 
least as fast as l/a 2 . 

B. Stationary state 

We now turn to systems already in the stationary state 
and give an overview of their statistical behaviour. 

The most striking feature of the stationary state are 
the large synchronized blocks. Sites within the blocks 
topple the same number of times, while sites closer to the 
boundaries topple less often. From time to time, a large 
avalanche that begins outside the block runs through 
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the entire block. Between the large avalanches, the sites 
within a block topple mostly one by one, lifting each other 
almost exactly to the threshold. Let us consider a small 
region within such a synchronized block and let us show 
that the sites in this region must have approximately the 
same height, given the dynamics just described. When 
the sites topple one by one, their height differences are 
exactly the same as before, after each site has toppled 
once. When an avalanche enters the region from outside 
and extends several sites beyond it, the change in height 
differences is calculated by multiplying the state vector x 
with the appropriate product S av of M* matrices. If the 
avalanche passes our region from the right to the left, the 
elements of S a " in a row i belonging to our region are 







for i 
else 



1 < 3 < 3i 



(4) 



with ji n i + 1 being the site that triggered the avalanche. 
If we denote with Xi and x\ the values of the force differ- 
ences within our region before and after the avalanche, 
we have 

x i = a(x' i+1 + Xi-i) < 2a Max (x' i+1 ,Xi-i) 

for an avalanche passing through the region from the 
right to the left. The asymptotic values of the xi after 
many avalanches satisfy Xi = a(xi + \ + within the 

synchronized region. This condition can only be satisfied 
with all Xi being zero or with Xi decreasing by a factor 
of the order a from one value of i to the next. Deep 
inside the synchronized region, the Xi become therefore 
very small. 




10 15 20 
site i 

FIG. 7: Thick line: mean number of topplings (t), as a func- 
tion of site index i, averaged over 1000 different synchronized 
systems and 10 4 topplings for a — 0.2. dotted line: analytical 
result for a = 0.2. g was choosen to be 0.6168576. 

Directly related to these blocks is the behaviour of (t)i, 
the mean number of topplings per site per unit time, 
which we observed as a function of the site index i, aver- 
aged over a long time and over many systems, as shown 
in figure [7| 



The fact that sites in the synchronized regions topple 
the same number of times, while those at the boundaries 
topple less often (due to the missing neighbours, or due 
to neighbours toppling less often), can also be explained 
analytically in the following way: At every site, a local 
balance equation has to be fulfilled. Let g be the mean 
force increment per unit time, and a the mean package 
size, which we assume to be the same at each site. The 
balance equation then is 
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for each site i, 

(*)0 = {t)L+l = 0. 

trix form 



with 



(t)i 



(t)i 



(5) 



Eq. 



the boundary condition 
JSJ can be written in ma- 



where T is tridiagonal and given by 



(6) 



LxL 



( 1 



V o 



-a 
1 -a 



-a 1 



1 / 



(7) 



LxL 



and dh is a vector with constant entries 1. The numer- 
ical solution of Eq. (J5J is also shown in Figure [7| An 
analytical solution can be found by making a continuum 
approximation to Eq. Ip|). 



d 2 {t{x)) 
dx 2 



1 



- 2 " Ul u 9 
- — ii - - 

a a 



(8) 



The solution that satisfies the boundary conditions is 



(t(x)) 



1 - 2a 
sinh (7(.x — L - 



(9) 



1)) 



sinh jx 



sinh( 7 (L + l)) sinh( 7 (i + l)) 



(with 7 = 



1-2S \ 



This mean-field result predicts that 



the thickness of the boundary layer is proportional to a, 
which agrees with our previous numerical finding that the 
time until the onset of synchronization is proportional to 
a. 

We also considered (s)i, the mean size of all the 
avalanches triggered at site i. The result for a = 0.2 
and L = 100, averaged over 10000 synchronized systems, 
is shown in figure [SJ Almost all of the large avalanches 
are triggered near the boundaries. Also shown in figure 
IHlis the relative number of avalanches triggered at site i, 
which also shows narrow peaks at the boundaries of the 
system, but also a broad peak in the center. 

Combining the two data sets, we arrive at the following 
scenario: In the stationary state, most of the avalanches 
are single topplings. All large avalanches are triggered 
near the boundaries and extend far into the synchro- 
nized block. If they do not reach the end of the syn- 
chronized block, the rest of the block topples in a se- 
ries of smaller avalanches, mostly of size 1. These small 



9 



30 rr 




FIG. 8: Solid line: (s)j: mean size of avalanches triggered 
at site i; dashed line: (# a }i: mean number of avalanches 
triggered at site i. Both data sets are plotted as function of 
the site index i, averaged over 10000 different synchronized 
systems for a = 0.2, L = 100. 

avalanches cause the broad peak at the center of Figure|Hl 
The structure of the peaks at the boundary of the curves 
in figure [3] depends on a, and results from the averaging 
over many different stationary states. 
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FIG. 9: Number of avalanches n(s) of size s, divided by 
the total number of topplings and the system size, averaged 
over at least 2000 systems Top left: L = 1000, a = 0.2. Top 
right: L — 500, a = 0.15; solid line: precision 2 -12 ; dashed 
line: precision 2 -20 ; dotted line: precision 2 -28 . Bottom left: 
L — 1000; solid line: a = 0.1; dashed line: a = 0.2; dotted 
line: a = 0.3; dash-dotted line: a = 0.4. Bottom right: 
a = 0.15; solid line: L = 100; dashed line: L — 500; dotted 
line: L = 1000. 

We now turn to the size distribution of the avalanches. 
For two-dimensional systems this is believed to obey a 
power law (see for example 0, H). However, the one- 
dimensional systems show this feature only for short 



times or for small system sizes. In the stationary state, 
the avalanche-size distribution looks like the upper left 
curve of figure^] It was obtained by averaging over 2140 
different systems and 10 9 topplings, where we neglected 
the first 10 9 transient topplings of each system. 

For small avalanche sizes, a power law is visible, but 
only for a single decade and up to the sharp cutoff. We 
see a large gap, followed by peaks centered at system size 
and half the system size. The shape of n(s) for small s 
depends only on the value of a and on the precision used 
in the simulations, but not on L. We checked this by in- 
serting more sites into the synchronized blocks and com- 
paring the resulting avalanche-size distributions, which 
only differed for the large avalanches. (See also the bot- 
tom right curve of figure |5J) . 

For larger precision, more avalanches are found. (See 
the upper right curves of figure |5J. The reason for this is 
that the period of the stationary state is longer, as stated 
already before. Smaller values of a result also in smaller 
avalanches, because a higher precision is needed in order 
to resolve force differences (which scale with powers of 
a). Note also the nonvanishing weight for avalanches of 
size 2L for a > a c in the lower left figure. 

The results for n(s) confirm the picture that the sys- 
tem is composed of a boundary layer that controls the 
dynamics and determines the stationary state, and a syn- 
chronized block of sites that topple the same number of 
times and that can be made larger without modifying the 
boundaries. 



VI. DISCUSSION 

The investigation of the one-dimensional version of 
the self-organized critical earthquake models has revealed 
many intreaguing features. When viewed as a dynamical 
system, the system shows 4 different types of attractors, 
all of them being periodic. In contrast to a 2-site ver- 
sion of the model, where the variables can only change 
continuously in time [TT| . and to a many-site version, 
where the reset rule is — > z% — z c |12| . we do not find 
chaotic attractors. In contrast to the one-dimensional 
Zhang-model, which is conservative and has a stochastic 
force input jl3 , the phase space volume does not neces- 
sarily shrink for systems that have the same sequence of 
topplings and avalanches. 

In the stationary state, the model consists of two 
boundary layers, the thickness of which is larger for 
longer attractors, and an inner part consisting of one 
or two synchronized blocks, where all sites have ap- 
proximately the same force value. The synchronized 
blocks can be made larger without changing the dynam- 
ics of the boundary layer or the period of the attractor. 
Large avalanches are always triggered near the boundary. 
These features are clearly reflected in the avalanche-size 
distribution, where the small avalanches are independent 
of the system size for sufficiently large systems, while the 
large avalanches are proportional to it. 
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Several features of the one-dimensional model are very 
similar to the two-dimensional model, while others dif- 
fer. In both versions, large avalanches are only triggered 
at the boundaries |5|, and synchronization proceeds in- 
wards according to a power-law in time |2lll4|. The inner 
part is dominated by avalanches of size 1 jJ0>^3] even 
in the stationary state. The computing precision affects 
the avalanche-size distribution ||. However, while there 
are less large avalanches for smaller computing precision 
in one dimension, there are more large avalanches in two 
dimensions. That the inner part can be made larger with- 
out changing the dynamics of the boundary region, must 
be a special property of the one-dimensional system due 
to the fact that the boundary of a synchronized block is 
merely a point and that avalanches can propagate only 
along lines. Nevertheless, the inner part is to some ex- 
tent slaved to the boundary region also in two dimen- 
sions. The exact interplay between the two still needs to 



be clarified. However, we want to conjecture that in two 
dimensions the avalanche-size distribution separates also 
in two parts, when the system size is only made large 
enough. We expect the first part to become essentially 
independent of the system size, while the second part be- 
comes proportional to it. However, in order to verify (or 
falsify) this conjecture, larger and faster computer simu- 
lations are needed than those that have been performed 
up to now. 

Thus, it appears that the OFC-earthquake model 
is not self-organized critical in the sense of exhibiting 
avalanches of all sizes. Rather, the avalanche-size distri- 
bution results from the combined effect of several mech- 
anisms, and only for sufficiently large system sizes do 
different types of avalanches become clearly separated in 
size. 

This work was supported by the Deutsche Forschungs- 
gemeinschaft (DFG) under the contract (Dr300/3-l). 
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